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I. INTRODUCTION 

Progress in laser technologies has resulted in the op- 
portunity to create ultra-strong electromagnetic fields 
in tightly focused laser beams. In the present paper 
we discuss the numerical methods designed to simulate 
processes in strong pulsed laser fields interacting with 
plasma. Attention is paid to the recently achieved range 
of intensities, J > 2 • lO^^W/cm^ [l[, and the larger in- 
tensities projected, J ^ lO^^W/cm^ [Ij. 

For a typical laser wavelength, A ~ l/xm, electron mo- 
tion in laser fields at J > lO^^W/cm^ is relativistic: 



|a| » 1, 



eA 



(1) 



where and e = — |e| are the mass and the electric 
charge of electron. 

However, if we want to evaluate the properties of an 
electron in a strong field as an emitting particle, more- 
over, a particle, which emits photons we need to be guided 
by the more severe condition: 



alal > 1, 



of energy, 2Ry = a^TOgC^. For a real laser, in addition to 
the field magnitude, importance rests on the laser photon 
energy normalized by 2Ry: 



huj 45nm 
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27.3eV. 



(3) 

For the majority of ultra-strong lasers this parameter is 
of the order of 10~^: W2Ry ^ 0.04 for the Nd-glass laser 
(A ~ 1.06 /im), W2Ry ~ 0.06 for the Ti-sapphire laser 
(A « 0.8 /im). Therefore, the following product. 



W2Ry 



da 



cEl 
47r 



Jp = « 2.4 • lO^^W/cm^ 



(4) 

is less than one even for planned intensities (although 
\da./d^\a might be greater than one). Herewith, the 
estimates are made for the ID wave field, a = a(^), 
^ = u{t—x/c), < ^ < ^max- Eq.(|31) is expressed in terms 
of the local instantaneous intensity of the laser wave, J. 
Note that the Left Hand Side (LHS) of Eq.d!]) equals the 
ratio, jEj/i^p, of the wave electric field, |E| = \/47r J/c, 
to the characteristic field, Ep = |e|/A^, constructed from 



(2) an elementary charge and the Compton wavelength: 



in which the fundamental fine structure constant is 
present, a = /{ch) ~ 1/137, linking its radiation to 
its motion (herewith, the subscript a denotes the di- 
mensionless parameter multiplied by a). With the re- 
cently achieved intensity of J '-^ 2 • 10^^ W/cm ~ 
(1/a^) • 10^^ W/cm^, this newly important dimension- 
less parameter exceeded unity! 

However, this estimate could be applicable only to 
a 'theoretical laser', for which the photon energy, fiw, 
would be comparable to the 'characteristic' atomic unit 



Ac = w 3.9-10-"cm, 



Ar 



27rAc ~ 2.4-10"^°cm. 



This field strength is associated with the Coulomb field 
between the components of a virtual electron-positron 
pair (which are 'separated' by the Compton wavelength) . 
Across the interval of intensities bounded by Ineq.dD) 
from below and by Eq.(|3|) from above, that is, at 



45nm 



W2Ry 




(5) 
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the role of important physical effects changes dramati- 
cally, incorporating radiation and its back-reaction, and 
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QED effects of electron recoil and spin as well as pair 
production. Given that currently available laser intensi- 
ties can access this kind of interaction, it is clear that the 
development of a suitable model is timely. 

Radiation-dominated laser fields. An accelerated 
electron in a strong laser field emits high-frequency radi- 
ation. The radiation back-reaction decelerates such an 
electron, the effect being more pronounced for longer 
laser pulses Q. In Ref.Q a condition for the field to 
be radiation-dominated is formulated in terms of the 
ratio between the magnitudes of the Lorentz force and of 
the radiation force, which gives: 



da 
d^ 



> 1. 



(6) 



Herewith the electron dimensionless energy, £, and its 
momentum, p, are related to mgC^ , and nieC correspond- 
ingly, and subscript |j herewith denotes the vector pro- 
jection on the direction of the wave propagation. 

While a strong laser pulse interacts with energetic elec- 
trons, which move oppositely to the direction of the pulse 
propagation, the condition £ — sa 2f ^ 1, facilitates 
the fulfillment of Ineq.®- In the course of a strong laser 
pulse interacting with a dense plasma the counterprop- 
agating electrons may be accelerated in a backward di- 
rection by a charge separation field. For this reason, the 
radiation effects in the course of laser-plasma interaction 
are widely investigated (see Refs.[3-Q) and efficient com- 
putational tools are in demand. 

QED-strong laser fields. In Quantum Electro- 
Dynamics (QED) an electric field should be treated as 
strong if it exceeds the Schwinger limit: |E| > Es = 
meC^ / {\e\Xc) (see Ref.Q)- Such field is potentially ca- 
pable of separating a virtual electron-positron pair pro- 
viding an energy, which exceeds the electron rest mass 
energy, mf,c^, to a charge, e = over an acceleration 
length as small as the Compton wavelength. Typical ef- 
fects in QED-strong fields are high-energy photon emis- 
sion from electrons or positrons and electron-positron 
pair creation from high-energy photons (see Refs . H-flOj ) . 

Here we assume that the field invariants (see [IJ]) are 
small as compared to the Schwinger field: 



(E.B)|«i?2, 



(7) 



B being the magnetic field. Below, the term "QED- 
strong field" is only applied to the field experienced by 
a particle. For example, a particle in the ID wave field, 
may experience a QED-strong field, Eq = \dA/d£^\uj{£ — 
P\\)/c, because the laser frequency is Doppler upshifted 
in the comoving frame of reference. The Lorentz- 
transformcd field exceeds the Schwinger limit, if 



X=lEo/Es = ^i£-p\0 



da 



» 1, 



(8) 



where X = c/uj. Numerically, the parameter, x, equals: 



da 
d^ 



0.7 



(g-P||) 

103 



J 



1023 [W/C 



Estimates for laser-driven electrons. In the crit- 
ical parameters as in Ineas. (|6l8p the factor, £ — is 
not linked to the wave intensity in the case where elec- 
trons are accelerated by an external source. In the course 
of the laser-plasma interaction, however, for bulk elec- 
trons (£ — p||) ~ £ ^ |p_l|- As long as the radiation 
back-reaction docs not dominate, the conservation law for 
the generalized momentum of electron gives: pj^ sa —a, 
and the LHS of Eqs.® may be evaluated in terms of 
the vector potential amplitude, Og. The wave becomes 
radiation-dominated (see IneqlH]), if: 



-1/3 



J > Jp(w2Ry)''/^ ^ (3-5)-1023w/cm^ 

Less certain is the estimate for the significance of QED 
effects. On one hand, for fields just approaching the 
radiation-dominated regime QED effects are already not 
fairly neglected: x ^ (3/2)(a;2Ry)^/3 w (0.5 - 0.6). On 
the other hand, in radiation-dominated fields the esti- 
mate for £ that we used above is not reliable. Because 
of this complexity, we surmise that the significance of 
QED effects in this regime can only be verified by direct 
numerical simulations. 

Paper content and structure. Numerical simu- 
lations of laser-plasma interactions become increasingly 
complicated while proceeding to higher intensities. At 
intensities J > 2 • lO^^W/cm the model should incorpo- 
rate the radiation back-reaction on the emitting electron. 
In this range x ^ 1 for bulk electrons, making the QED- 
effccts negligible. This model is presented in Sec. II. At 
J > 3-10^3\Y/cm^ the QED corrrections should be incor- 
porated to achieve a quantitative accuracy for electrons 
with X ^ 1- These corrections may be found in Sec. HI. 
At larger intensities, J > lO^^'W/cm^ the high-energy 
photons emitted by electrons and positrons produce a 
macroscopically large number of electron-positron pairs, 
as shown in Sec. IV. 

In each section we first summarize the theoretical 
model. Then we provide the analytical solutions, which 
may be used to benchmark numerical models. After this 
we describe the elements of the numerical scheme. 



II. QED-WEAK FIELDS (x < 1) 
A. Theoretical notes 

1. Emission spectrum 

In Ref.[ll| the spectral and angular distribution, 
d£T:iid/{di-o'dn), of the radiation energy, emitted by an 
electron and related to the interval of frequency, dio' , and 
to the element of solid angle, dn, for a polarization vec- 
tor, 1, is described with the following formula: 



rf£'rad(a^',n, 1) 
doj'dn 



i(Aci(a^')-i*)r 



(9) 
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Here the vector amplitude of emission, Ac\{u>'), is given 
by the following equation: 



Aci(a;',n) 



v(t) exp ( iJ \t - \ 1 dt. 



see Eq. (14.67) in Ref.[ll[. We express dS^;^^/ [dio' dn) 
in terms of the time integral of the radiation loss rate, 
dlci/{du}'dn), which is related to the unit of time, the 
element of a solid angle, and the frequency interval, and 
is summed over polarizations: 



d_ 

dt 



E 



d£y, 



ad 



dbj'dn 



dut'dn 



The spectral and angular distribution of the radiation 
loss rate is given by the Fourier integral: 



duj'dn 4:Tt'^c£{t) 



+ 00 



•p T - 



c 



X exp 



H 



-C/2 



[k' ■p{T')]dT' } d(. 



The cogent feature of the particle relativistic motion 
in strong laser fields is that the emitted radiation is 
abruptly beamed about the direction of the velocity vec- 
tor, p(T)/|p(r)|. Therefore, the angular spectrum of 
emission can be approximated with the Dirac function: 



dlciir) 



dui'dn \ IpI / dw' 

In the frequency spectrum of emission, 



dlcijr) 
doj' 



X sm 



P T 



■p\T 




T+C/2 



-c/2 



{[p{T)-p{T')]-l}dT' 



dc, 



for relativistically strong wave field, satisfying Eq.([T]), the 
sine function oscillates fast, so that the main contribution 
to the integral determining the emission spectrum comes 
from a brief time interval with small values of C, resulting 
in a universal emission spectrum: 



dlciir) ^ Qei(r) 



dio' 



Id 



2e^(/Le-/Le) 



(10) 



Here Qci(^) is the unity- normalized spectrum of the gy- 
rosynchrotron emission, such that j Qc\{''^)dr = 1 and 
K^{r) are MacDonald functions. We use the dimension- 
less photon frequency, w' = Tiuj' /(meC^), the character- 
istic frequency, Hic = fitUc/ (jrieC^), and the dimensionless 
wave vector, k' — Kk! , for emitted 7-photons and omit 
tildes in the formulae. Both the radiation energy loss 
rate, /d, and the QED-strength parameter. 



X = 7; 



3 y^C\/ —{Ihe ■ jhe) 

2 



(13) 



are expressed in terms of the 4-square of the Lorcntz 4- 
force: /j^^ = Ei^Le. ■ u/c, fie), where u = cp/ v^l + 
is the velocity and f^e = eE + |[u x B] is the Lorentz 
three-force. 

Thus, the acceleration of electrons by the laser pulse 
must be accompanied by gyrosynchrotron-like emission 
spectrum (which is actually observed - see Refs.fl^ [isj). 
The general character of such emission spectrum had 
been noted in Re f.|12| (this comment may be also found 
in Sec. 77 in Ref . [14| ) . The material of the present sub- 
section was published in Ref. [isj . 



2. Equation for the radiation emission and transport 

The above considerations justify the method for calcu- 
lating the high frequency emission as described in Ref. Q 
(see also [13] )■ In addition to calculating the electro- 
magnetic fields on the grid using a Particle-In-Cell (PIC) 
scheme, one can also account for the higher-frequency 
(subgrid) emission, by calculating its instantaneously ra- 
diated spectrum. Compared with direct calculation of 
the RHS of Eq.® (the means of calculating the emis- 
sion used, for example, in Ref. [lH . [l7| ) the approach sug- 
gested here, though mathematically equivalent, may be 
decidedly more efficient. 

Generally, the Radiation Transport Equation (RTE, 
cf. Ref.ll^) should be solved for the radiation energy 
density, related to the element of volume, dV (herewith 
the symbol is omitted): 



d£, 



rad 



duj'dndV 



An electron located at the point, x = Xe(i), contributes 
to the right hand side (RHS) of the RTE as follows: 



dl 
'dt 



c(n • V)I 



^ciQci(^) ^2 



Qy/S f°° / / oj' The LHS of the RTE describes the radiation transport, 

~ "StT^ / ^^/^^^ ^'^^ ' ^ ~ ^^'^^ while in the RHS, in addition to the emission source, 

^ there should be the terms accounting for the radiation 

absorption and scattering. However, at x ^ 1 and at 
Wc = Sx- (12) realistic plasma densities these effects may be neglected. 
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Under these circumstances the RTE can be easily inte- 
grated over time and space, giving: 



d8. 



ad 



duj'dn 




Wcl(r)j2/„ P 



P / 



(14) 



Since Ea. (jl4p depends on frequency only via Qcii'f), one 
can calculate instead of Eq. an integral as follows 



d£. 



(m) 
rad 



dndcu 



f \y.^-^6^ fn-P)5(logc^-log 



UJc 



dt. 



Once this modified spectrum has been integrated over 
the whole simulation time, a true spectral distribution 
can be recovered using a simple convolution as follows: 

dnduj' J \Lij J dndu) 

Manifestly, the result is the same, which allows one to 
avoid calculating the spectrum, Qdir), at each time step. 



All 4- vector equations in this paragraph are written with- 
out indices and the tensor multiplication is denoted with 
dot-product and/or powers of tensor, e.g. F -u ~ F'^'^Uk, 
F'^ ■ u = F ■ F ■ u ~ F^^Fmu^ etc. Now we re-write the 
LL equation: 

du \ ~ / ~ \ dF u , ~, ^ , . 

— ^ —F ■ [u + F -u] + — ■ u (u-F^-u). (19 

dr tq \ I dr tq 

and compare it with Eas. (|16ll7p : 

'^^ = Lpu- —{p- F'^ -p), u=p + F-p. (20) 
dr To To 

Solving the momentum from the second of Eqs.ipUl). 
P ~ X]^o (~^)" ■ Accounting for the anti-symmetry 



of the field tensor, the first of Eqs.([20|) may be re- written 
for 4- velocity and 4-acceleration, similarly to Eq. ([T9]) : 



du 1 ~ , „ 
— = —F- [u + F -u 
dr To 



dF 



^.^(-F)"..-^(..^P"..).(21) 



To 



3. Equation for electron motion: accounting the radiation 
hack-reaction 



Here we use the equation of motion for a radiating 
electron as derived in Refs.j^, [l^. In 4- vector form this 
equation may be written for the electron 4-momentum, 
p", normalized per mgC, in terms of the Lorentz 4-force, 
= eF^Ppp and the field 4-tcnsor, F"'^ = (E,B): 



dp"' 



:F' 



dr 



/clP° 



dx" 
~d^ 



cp 



2e' 



To 



Three- vector formulation of Eas. (|16lir7)) is: 



dp 

'di 



he _^ e[u X B] ^ Ug^(u ■ fLe) 



(17) 



(18) 



It 



u + u, u 



To f^e ~ U(U ■ fLe)/c^ 

nie l + To{u-iLe)/{mc^)' 



4- Comparison with the Landau- Lifshitz equation. 

Many authors simulate motion of an emitting electron 
using the equation as suggested by Landau and Lifshitz 
(LL - see Eq.(76.3) in [3]); motivating a comparison with 
the approach we use. To simplify the formulae, we intro- 
duce the 4- velocity, and normalize the field tensor: 



1 dx" 
c dr 



ToeF 



'ik 



2aF 



ik 



3Es 



The only distinction from Eq.(|T9]) is that in Eq.([2T|) the 
infinite series are present, while in Ea.([T^ there are only 
starting terms of these series. 

How large is the difference numerically^ For the second 
series one may evaluate both the total sum: 



oo 

E 

n=l 



F'^'^-u 



p- F -p^ \-a 



9^ 



(16) and the residual sum, omitted in Eq.([T 



n=2 



9-^ £;? 



E^ 



The residual sum is reduced by a factor, {2a/i)'^{E'^ — 
B'^)/Eg, which is small according to Ineq.(l7]). 

How theoretically important is the distinction be- 
tween the two approaches? We discussed this issue in 
Refs.jl, [l^ and noted that the LL equation conserves 
neither the generalized momentum of electron nor the 
total energy-momentum of the system consisting of an 
emitting electron, the external field and the radiation. 
Another distinction is that the LL approach maintains 
the identity, = 1, while Eqs. (|20[) maintain more im- 
portant identity, = 1. turning to the Dirac equation 
in the limit of QED-strong field. For the square of the 
4- velocity we obtain: 



p-F^ -pKil- 1.05 • 10- 



(22) 



which is not exactly unity, but in QED-weak fields, x 
1, the distinction is negligible. 

The computational advanta,gc of Eas. (|20|) as compared 
to the LL equation, is, first, the higher efficiency: com- 
pare the compact expression for three- force in Eq. (|18p 
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with that given in [iJl (see section 77, problem 2). Sec- 
ond, the numerical scheme for Eg. ljlSp is more reliable, as 
it is bound to yield total energy conservation. Thus, even 
for fields in the QED weak regime, the use of Eqs. ipD)) is 
more suitable than the use of the LL equation. 



B. Analytical solution 

Pertaining to the validation against a semi-analytical 
theory, we begin by describing the spectrum of emission 
from an electron in the field of a ID circularly polarized 
wave. A constant wave amplitude, Qq is assumed to be 
below the radiation-dominated regime. In this case p_L ~ 
—a, so that f ^ — = 1 -|- Og. The modified spectrum 
can be expressed in terms of the characteristic frequency, 
which is a function of the current value of the electron 
energy: 



UJcO 



1 



(g-P||)^ 

1 + 



(23) 



as well as the parameter, ujco which is introduced as the 
following function of the wave amplitude and frequency: 



UJcO = -^2YiyC? (ao -I- aj^) 



(24) 



The maximum frequency of emission is determined by 
the initial momentum of electron: 



1 



Then, ^*($) is a normalized phase: 



^0 



(25) 



For the whole pulse the total normalized phase 



Coo = g^^aoy^l + alUe.^ -^yj 



2 / J \ f 



(w2Ry)2 ' 



characterizes the capability of the pulse to arrest the 
counterpropagating electron by means of the radiation 
back-reaction. Particularly, a pulse of duration corre- 
sponding to > 1 arrests an electron of any energy. 
The modified spectrum has a shape close to a power-law 
(see derivation details in [20j): 



dE, 



(rn) 
rad 



doj 



4w, 



cO (w/wco-l) 



3/2 ' 



(26) 



^cmin <^ ^ <^ ^cmax 



WcO (^cQ ^cO 

where Wcmin should be found from Eg. ip^ . for given <^^: 

/ \ -2 



1 



(28) 




1 10 100 

Frequency 

FIG. 1; The shape of normalized spectra, [dSrad/ dui') ■ 
[^LOco I {irieC^y/l -I- a^)], versus the normalized frequency, 
ui'/ujco, for different pulse duration. The figure is scalable, 
particular choice of physical parameters, may be the follow- 
ing: |a| = 50, £ = 180 MeV, pulse durations are: 5 fs (curve 
1), 36 fs (curve 2) and 220 fs (curve 3). The spectrum broad- 
ening and softenning is due to the radiation reaction. In the 
absence of this reaction curve 1 without changing its shape 
would scale proportionally to the pulse duration. A zero value 
of logcj' corresponds to ~ 150 keV. 



The true (transfromed) spectrum can be obtained from 
the modified spectrum as in Eq. (|26l) by applying a convo- 
lution transformation following Ea. (|15p . The longer the 
pulse, the more softened and broadened the radiation 
spectrum becomes (see Fig.l). 



C. Numerical model 

Now we introduce the following normalized variables: 

u 

c 



E 



t — LUt, 

|e|E ^ 

meCUJ ' 



X ~ wx/c, u 



B 



\e\B ^ 

m.f.cLO ' 



47r|e|j 



Note that the electric current density, j is normalized 
per |e|ncrC, where Ucr = mea;^/(47re^) is the critical den- 
sity. Below, we use these dimensionless variables and 
omit tildes in notations. The equations of motion for 
electrons and positrons read: 

~ ^Le.p "F [Ue,p ^ B] — Up p£^ p (fie.p ' Ue.p) j 



dt 



Pip, 



Pe,p 



(27) normalized Lorentz force and u being: 



fie.p = T (E + [Ue_p X B]) 



TqUJ 



1 + ^(Ue^p ■ fLe,p) 



(29) 
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For reference we also provide the energy equation: 



= T ([Ue,p + U,,p] • E) - Ei^ (fLe,p ' Ue,p) , (30) 

For ions with the charge number, Z, and the mass, Afi, 
the momentum is normahzed per M^c, so that in their 
equation of motion the electron-to-ion mass ratio comes: 



dp, Z 



dt Mt/m., 



(E+[u, X B]) 



(31) 



Eas. (130l33l34D give: d{£pi^ 



£ra.d) /dt = 0, where 



d£, 



rad 



dt 



(35) 



is the radiation energy loss rate. Therefore, the contri- 
bution from electrons and positrons to the dimensionless 
radiation energy at each time step. At, is calculated as 
^e.p (fie,p • Ue,p) At/Nppc- OncB integrated over the sim- 
ulation time, the radiation energy may be converted to 
physical units on multplying it by a factor of Sq. 



dX4 

dt 



U,; 



(32) 



Below we assume that Z = 1 and Mi = Mp is the proton 
mass. The normalized Maxwell equations read: 



5E 
'dt 



V X B, 



SB 

'dt 



-V X E. 



(33) 



Macroparticles and their current 



We assume that the rectangular grid splits the compu- 
tational domain into the control volumes (cells), AV^ = 
J^Axfc. If the electron density equals ricr, there are 
UcrAV electrons per cell. The latter number is typically 
very large, so that the plasma electrons cannot be simu- 
lated individually and they are combined into macropar- 
ticles with a large number of "electrons-pcr-particle" , 
iVepp. In a plasma of critical density the number of 
(macro)particles per cell is: Npp^ = nci-AV/Nf,pp- 

The electron current density inside the given cell is ex- 
pressed in terms of the sum over electron macroparticles 
in this cell. As long as the electric current density is nor- 
malized per ricrC, the contribution to the latter sum from 
each macroparticle equals —{l/Nppc)dx.e/dt. On adding 
the contributions from positrons and protons, we obtain: 



Ee (Ue + Ue) + Ep ("p + "p) + E, 



(34) 



2. Energy integral and energy balance. 

We now establish the relationship between the energy 
integral and the finite sum which represents this inte- 
grals in simulations. Particularly, the field energy may 
be calculated as the total of point-wise field magnitudes 
squared: £fioid = ^Ece/is (^^ +B^), which is by a fac- 
tor of £o = WeC^ (ricr Ay) different from the dimensional 
field energy. Now consider the total plasma energy, which 
includes the particle energy as well: 



EceHs ( ml Ei + Ee + Ep 



-'plasma 



-field 



N, 



3. Algorithmic implementation 

The algorithmic changes to the standard PIC scheme 
are minimal as long as we ignore the radiation trans- 
port and only integrate over time the energy emit- 
ted by electrons (and positrons, if any). To collect 
the radiation, we introduce the energy bins (an ar- 
ray) E-cadijk{\og{uj)i,9j,ipk); which discretize the modi- 
fied frequency-angular spectrum of emission. Inside the 
desired interval of the photon energies we introduce the 
logarithmic grid, log(a')i, equally spaced with the step, 
A log(w). We also introduce a grid, Oj, (pk, for the two po- 
lar angles of the spherical coordinate system, with Arij^ 
being the element of solid angle: An^fc = sm{9j)A9Aip. 

To calculate both the spectrum of emission and the 
radiation back-reaction we modify only that part of 
the PIC algorithm which accounts for the electron mo- 
tion. Specifically, we employ the standard leapfrog nu- 
merical scheme which involves, among others, the fol- 
lowing stages: (1) for each electron macroparticle, up- 
date momentum through the time step by adding the 
Lorentz force, following the Boris scheme: p"^^''^ = 
p""^/^ + At fLe(E",B"); (2) solve the energy and 
the velocity from the updated momentum: = 
Vl + (p"+i/2)2, ur^/' = p:+'/^ (3) use the 
calculated velocity to update the particle coordinates: 
x"+^ = x" + Ue^^^^At and account for the contribu- 
tion of the electric current element, —u'e'^^^^ At, to the 
Maxwell equations. Again, these stages are standard and 
may be found in [2l| . We introduced new steps into this 
algorithm between stages (2) and (3) as follows. 

2.1 Once stage (2) is done, recover the Lorentz force: 



(P 



n+l/2 



)/At. 



2.2 Find X = |w2Rya'^^"+^/VfL " i^Le ■ U"+l/2)2. 

2.3 Find Ue by putting f^e and Ue~^^^^ into Eq. ([29| . 

2.4 Calculate Wc = £e~^^^^X ^nd find the discrete value 
of log(ct')i most close to logWc. Find the angles, 
6j,ipk closest to the direction of p"+i/2. Add the 
radiation energy into the proper bin 



ppc 



(g"+V^)^(fLe-Ue)At 

WcA^ppcA log(w)Anjfc ' 
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FIG. 2: Test simulation result: (a): radiation energy 
spectrum (line 1), d£^s_d/duj' , and the modified spectrum, 
d£i!!,d/diu (line 2); (b),(c): the angular distribution of the ra- 
diation at instants: (b) t=lT, (c) t=100T, where T = 2tt/u}. 



n+l/2 
Pe 



n+l/2 
Pe 



2.5 Add the radiation force: 

At {-[Ue X B"] - ur'/'(fr+'/')2(fLe • Ue)}. 



n T^' 1 n+l/2 n+l/2 / /i , / n+ 

2.6 Find Ue ' = Pe ' /\Jl + {pe 

use this velocity through stage (3). 



-1/2 



+ Ue and 



Note that the algorithm modification is applied only to 
electrons (positrons), keeping unchanged the ion motion 
as well as the fields. 

The frequency-angular spectrum may be reduced to 



the frequency one: 8i- 
gular one: £, 



ad i 



.An,fc, to the an- 



rad jk 



I]j'?radiifct^iAlog(w) or to the total 



radiation energy: f^ad = Y^nk fradij7cW.jAlog(a;)An^ 



While postprocessing the results, we apply the con- 
volution transformation, ([TS|. to the radiation spec- 
trum and multiply it by Eq. The resulting spectrum, 
d£i-a,d/{dndu!'), is a function of log[7ia;'/(meC^)]. 



D. Simulation result 



The analytical solution presented in Sec lIIBI has been 
used to benchmark the numerical scheme. In the test 
simulation electrons with an initial momentum, = 300, 
propagate toward the laser pulse with sharp (2A) fronts. 
The circularly polarized laser pulse has amplitude, |a| = 
15, and duration, 100(27r/u;). Interacting with the pulse, 
the particles radiate energy, finally approaching momen- 
tum of p|[ w 130. 

In Figl^l^a) the spectrum of the resulting radiation 
dSi-nd/dio' is shown. We also provide the modified spec- 
trum (the distribution over uj = uic), which is close to 



satisfying a power law, and in a full agreement with the 
analytical solution. In FigI2][b)-(c) typical evolution of 
the angular radiation distribution, d£ra.d/dn, is provided 
for the same simulation. One can see that the majority 
of the radiation is concentrated in a narrow angle with 
respect to the direction of backscattered light (0°). A 
softer part of the radiation exhibits a wider angular dis- 
tribution and becomes less intense [Fig. [D^c)]. 



III. QED-MODERATE FIELDS (x ~ 1) 

When X is not small (J ~ 3 • lO^W/cm^), QED ef- 
fects come into play. Here we describe how to extend the 
methods used above towards finite x- This is achieved by 
applying realistic QED spectra of emission as derived in 
Ref.[15[, with the radiation force modified accordingly. 

This approach is applicable as long as we ignore the on- 
set of some new effects which are only pertinent to QED- 
strong fields. Specifically, while employing the radiation 
force, dp'^^^^/dr, it is admitted that the change in the elec- 



dp'' 



rad 



/dr, within the infinitesimal 



tron momentum, dr 
time interval, dr, is also infinitesimal. This 'Newton's 
law' approximation is pertinent to classical physics and 
it ignores the point that the change in the electron mo- 
mentum at X ^ 1 is essentially finite because of the fi- 
nite momentum of emitted photon. The approximation, 
however, is highly efficient and allows one to avoid time- 
consuming statistical simulations. Its error tends to zero 
as X — >■ 0, and it is sufficiently small at x 1- 

Another effect which we ignore in this section is the 
pair production due to 7-photon absorption in the strong 
laser field. This neglect allows us to avoid solving the 
computationally intense radiation transport problem. 



A. Theoretical notes 

1. Emission spectrum 

The emission probability found in Ref.flU within the 
framework of QED can be reformulated in a form similar 
to Eq.(l9]). The polarized part of emission may be reduced 
to Eq.Q with the modified vector amplitude: 



Aqed(w') = Jtt-AcI {-^ 



C 



{k-pfY 



subscript i and / denoting the parameters of an electron 
prior to and after the emission of a single photon, and: 



duj' 



C 



duj' 



Within the framework of QED the electron possesses not 
only an electric charge, but also a magnetic moment as- 
sociated with its spin. Usually the spin is assumed to be 
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FIG. 3: Emission spectra for various values of x- 

2 t 
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l°gio(4i/^c) 

FIG. 4: Emitted radiation power in the QED approach vs 
classical (solid); an interpolation formula /qed = /ci/(l + 
1.04^/ci//c)''''' (dashed). Here, Ic = /ci/x'- 



2. Equation for electron motion: accounting the radiation 
back-reaction. 



As long as the QED effects modify emission, 



Qciir) (3qed(7-,x), 



'QED, 



(38) 



the radiation back-reaction needs to be revised accord- 
ingly. In Refs.fsl. IT9| it was noted, that QED is not com- 
patible with the traditional approach to the radiation 
force in classical electrodynamics, while Eqs. (|16ll7|) may 
be employed at finite value of x on substituting /qed 
for Id- Alternatively, within the framework of QED the 
radiation back-reaction may be found by integrating the 
4-momentum carried away with the emitted photons and 
that absorbed from the external electromagnetic field in 
the course of emission. For a ID wave field this procedure 
gives the following equation (see Ref . [l5l| ) : 



dr 



QED 



TTleC TTleC 



(39) 



and in such field F"^^ Fpf,p^' / {p^F'^P Fp^,p^') = k^/ik-p), 
k" being the wave 4- vector. Eg. ([59]) coincides with 
Eas. (|16ll7p . in which the substitution (p8| is made, or, 
the same, the three- vector formulation for u is applicable 
to the ID wave field, if the substitution is done as follows: 



To To 



/qed 



Toqix)- 



(40) 



depolarized (as is done in Ref. |T5| ) . and, accordingly, a 
depolarized contribution to the emission appears: 



dio' LJc Stt Cfi \Cfi 



Thus, the QED effect in the emission from an electron 
in a strong electromagnetic field reduces to a downshift 
in frequency accompanied by an extra contribution from 
the magnetic moment of electron. The universal emission 
spectrum in QED-strong fields is given by the total of 
Eqs.dSSET]): 

^^^TT^ — 9(x)<3qed(?',x), /qed = IcMx), 

where Qqed = Qqed/? the normalized by unity spec- 
trum, q{x) = /o°° Qqed(^' x)'^''' ^ 1 is the normalization 
parameter, the spectrum before normalization is: 



<3qed(^,x) = ^z—"^ 



Stt 



and — r/(l — x^). 



K5/3{r')dr' + x'rr^Kysirx) 



Although this approach is derived for the ID field, we 
may apply it to an arbitrary 3D focused field. An argu- 
ment in favor of this generalization is that the property 
of a ID wave, {k ■ k) = 0, which is used while deriving 
Ea. (p9)) . holds as an approximation for any field. Indeed, 
on calculating the 4-square of F°'^Fp^p'^/ {p^F'^'^ Fpf^.p'^), 
which 4-square is similar to (k ■ k)/(k ■ p)'^ , we find: 



P 



pi 



■P 



9 



ip-F^-pY 4^ 



El 



81 

16' 



(E_B)2 



« 1, 



the inequality holds at x > 1 according to Ineq.^. 

Another criterion which should be checked at x > 1 is 
the requirement for the difference, (1 — m^) cx x^? ^-s in 



Eq.(I221) to be small, 
we find that (1- 



Applying the substitution ([iU)) to 
u^) < 2 • 10~^, reaching its maximal 
value at X ~ 3.4. This 'error' is negligible, even if one 
assumes, than any theory allowing 7^ 1, is erroneous. 



B. Analytical result 

In FigOwe show the emission spectrum for an electron 
interacting with a laser pulse (see [l5| for detail). We 
see that the QED effects essentially modify the spectrum 
even with laser intensities which are already achieved. 



-3-2-10 1 2 3 4 
l"gio(f-/Wco) 



IkeV 1 MeV 1 GeV 



FIG. 5: The emission spectrum for 600 MeV electrons inter- 
acting with 30- fs laser pulses of intensity 2 • 10^^ W/crn^: with 
(solid) or without (dashed) accounting for the QED effects. 
Here huco ~ 1.1 MeV for A = O.Sfj,m. 



C. Numerical model. 



FIG. 6: Backscattered light in simulation of the interaction of 
laser pulse of intensity 2 x lO^^W/cm^ with plasma of density 
6.5 X lO^'^cm"^. For each xi (vertical axis) the convolution 
integral (the term in the above formula for convolution) is 
calculated and presented as a function of htj' (horizontal axis). 
The line shows the total emitted energy as a function of the 
cutoff photon energy (i.e. the integral spectrum). 



As long as the QED spectrum of emission depends on 
X, the bins for collecting the radiation energy should be 
refined: £i-a.d ijkiii^i,'n.jk,xi)- Once for a given electron 
(or positron) the parameter x is calculated; the discrete 
value of xi should be found most close to x- Then, pa- 
rameter e should be found following Ea.(|40l). e = qiT^uj, 
using pre-tabulated value, qi = q{xi)- Then, Ug should 
be expressed in terms of e and the radiated energy should 
be added to a proper energy bin: 



ad ijkl 



^radijkl H~ 



uJcNppcA\og{Cj)An 



-jk 



While postprocessing the results, a convolution similar 
to (fTSj) should be applied with x;-dependent spectra: 



dgi.ad(^',n) 
dndu' 



fo^y Qqed {^,Xi^ ^^adijkidlogL 



D. Simulation results 

In a ID simulation presented in Fig|6]a linearly polar- 
ized laser pulse with a step-like profile having 2-A front 
and amplitude a = 300 interacts with plasma of density 
no = 30??c during 50 cycles. About half of laser energy is 
converted to high energy photons. The data for backscat- 
tered photons indicate that values of x ~ 1 are achieved. 
These values are reasonable, as the energy of electrons 
moving toward the pulse is as high as ISOMeV, and the 
vector potential approach values a = 400 due to super- 
position of incident and reflected light. One can see that 
~65% of emitted photons exceed lOMeV, and 96% are 
above IMeV. 



IV. QED-STRONG FIELDS (x > 1) 
A. Theoretical notes 

In cases where x ^ 1 one needs to solve the RTE 
in order to account for 7-photon absorption in strong 
fields. This may be done using the Monte-Carlo method, 
in which the radiation field is evaluated statistically (see 
Refs. [23 - [24| ). Instead of the radiation energy density 
the photon distribution function may be introduced as 
follows: 



/^(x,a; ,n) = 



(41) 



Similar to the way the electron macroparticles represent 
the electron distribution function, the photon marcopar- 
ticles may be employed to simulate the photon distri- 
bution function. To simulate emission, the photons are 
created with their momentum selected statistically. The 
photon propagation in the direction of n is simulated in 
the same way as for electrons. The absorption with the 
known probability is also simulated statistically. 

Now we may split the radiation for the part which may 
be treated in the way we followed so far (see Sec.II-III) 
and for the Monte- Carlo photons. We choose a param- 
eter, X* (0.05 — 1) and assume that: (1) an electon 
with X ^ X* contributes only to d£rad/{duj'dn); while 
(2) for an electron with x > X* the regular spectrum 
of emission, Q'qed (''i x) (which is normally truncated at 
r = I/x), is now truncated at r = X*/x^ < 1/x a-i^d 
the emission of photons with x*/x^ ^ ^ Xj O'^: the 
same, ^ < ^ < 1, is treated statistically. The regu- 
lar radiation loss rate as well as the contribution to the 
radiation force cx — p/qed should be both reduced by a 
factor of qt(x)/'z(x) a-t x > X*j where the truncated spec- 
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* / 2 

trum integral equals: qt{x) = Jq Qqed(''' x)^'"- The 
normalized by unity trucnated spectrum is: Qt{r,x) — 
'^'^^^Q™''''^^ r < x*/x^- The emission probability to be 
used at X > X*; ^ > X*/x^ niay be found in [25| : 



e— ^e,7 



where 



field we use an equation 



and for a ID wave 
If the 



V3ad(, 
2X 



emission probability is averaged over time or over an en- 
semble, we return to the above spectrum of emission: 
u'ruec^ (dWe^c-r/idrdt)) = /qedQqed(?', x)j Here we 

apply the formula, —t^ — „ ^^ft™ i , which is also 

used in the numerical scheme. 



B. Semi-analytical solution 

In we demonstrated that as long as the distribution 
functions, /e,p,7, for electrons, positrons and photons in 
a ID wave field are integrated over the transversal com- 
ponents of momentum, their evolution is described by 
simple kinetic equations with the collision integrals. We 
solved these equation numerically. This choice of initial 
conditions corresponds to the 46.6 GeV SLAG electron 
beam and the laser intensity of J w 5 • lO^^W/cm^ for 
A = 0.8/xm, to be achieved soon. As long as the Monte- 
Carlo method is not used, the numerical results, such as 
the total pair production, may be used to benchmark the 
numerical scheme described here. 



C. Numerical model 

The modification of the numerical scheme as used in 
Sec. Ill is needed only for electrons and positrons with 
X > X* ■ The radiation energy added to the proper energy 
bin is corrected as follows: 



^rad ijkl ^ ^rad ijkl 



q{xi) uJcNppcAlog{Ld)Anjk 



and the same correction factor is applied to the second 
term in braces at algorithm stage 2.5. After stage 2.5 
a probable hard photon emission from the electron with 
X > X* is accounted, using the probability: 



(5 = 



^2£;(n+l/2) 



9V3 



The total probability of emission is given by a com- 
plete integral: We^e,-f = S /^./x"^ (t'^) ^ {t)- 
within the QED perturbation theory and within the 
Monte-Carlo scheme W, 



is assumed to be less than 
one. The probability of no emission equals 1 — We^e,-y ^ 
0. The partial probability, We^e,j{(^' < Wg), for the 
emission with the photon energy not exceeding the given 
value, Wg, is given by the incomplete probability integral: 



e— fe,7 



X* Ix 



X 



Therefore, for given 5 and x and for a randomly generated 
number, < rnd < 1 one can solve lo' /£ from an integral 
equation as follows (see detail in Appendix B): 



'/£ 



x' Ix 



w (z, x) dz = —— < 



X' Ix 



w{z,x)dz, (42) 



if the gambled value of rnd does not exceed We^e,-y- 
< rnd < We^e.y Otherwise (if We^^-y < rnd < 1) 
the extra emission does not occur. With calculated uj' /£, 
the emission is accounted for by creating a new photon 
macroparticle with the momentum, p^J^~^^^'^\u!' /£) and 
an electron recoil is accounted for by reducing the elec- 

(n+l/2) (n+1/2),^ , i i c-w 

tron momentum, Pe -)> Pe (1 - {'■^ l^))- 



1. Photon propagation and absorption 

The new element of the numerical scheme is the photon 
macroparticle, which simulates {ric-cAV) / Nppc real pho- 
tons. Its propagation with dimensionless velocity equal 
to n is treated in the same way as for electrons and ions. 

If the photon escapes the computational domain, its 
energy should be accounted for while calculating the total 
emission from plasma. For this purpose we introduce 
the energy bins, £^^*°*^^(log(a;')i, 6*^, (^t), such that the 
logarithmic equally spaced grid for the photon energy, 
log(tt'')i, and the polar angle grid coincide with those 
introduced above. The contribution from the escaping 
photon with total energy, ui' meC^ricilSV / Nppc, should be 
added to the bin with the closest log(a;')i, (/Sfe, with 
the macroparticle energy being converted to the spectral 
energy density by dividing by Aw' = w'Alog(w'): 



p,(tot) 
^rad ijk 



{tot) ^ meC^{ncrAV) 

NppcA\og{u') 



rad ijk 



The photon absorption with electron-positron-pair cre- 
ation is gambled in the same way as the emission (see 
details in Appendix B). Other absorption mechanisms 
may be also included. In postprocessing the simulation 
results, the softer 7-photon emission should be added to 
the total radiation spectrum: 



'^^rad 

dndu) 



11 



XI 



it \ -r^,Xi] £radijkid\ogu; 



Using the known relationships, 2 '^^^^'^ = —K^-i(r) — 
Kv+i{r), K-^{r) = K^{r), more integrals may be re- 
duced to the MacDonald functions, particularly: 



D. Simulation result 

Repeating the test simulation as described in Sec. III. D 
and applying the Monte-Carlo scheme at x > X* = 0-1: 
and without photon absorption, we notice only an in- 
crease in the fluctuations of the high-energy portion of 
the radiation spectrum. In this region, the photons are 
statistically underrepresented, the number of particles 
per A log(a;') being small. 



+ 00 



1 + 2^2 



3 fz' 



+ 00 



dz^ I K5/3{r )dr 



The advantage of the MacDonald functions is the fast 
convergence of their integral representation: 



Ki,{r) = / cxp [— r cosh(2)] cosh(t^z)(iz, 
Jo 



V. CONCLUSION 

Thus, the range of field intensities which may be sim- 
ulated with good accuracy using the described tools is 
now extended towards the intensities as high as (2 — 3) • 
lO^^W/cm". In such fields, which are typical for the 
radiation-dominated regime of the laser-plasma interac- 
tion, the suggested scheme is validated against a semi- 
analytical solution. Different versions of the equation 
of the emitting particle motion are compared and their 
proximity is demonstrated. 

Extension of the model for QED-strong fields can be 
easily incorporated into the scheme. The emission spec- 
tra are substantially modified by QED effects and simu- 
lation results for a realistic laser-plasma interaction are 
provided. 

For the QED-strong field regime of laser-plasma inter- 
action the Monte-Carlo method should be used to simu- 
late emission-absorption of harder 7-photons. Although 
such simulations are doable (see [26|), more work on the 
scheme validation is needed. 



Appendix A. MacDonald functions 

The MacDonald functions allow solutions for the fol- 
lowing integrals; 



+00 





■3 






/cos 


-r 
2 




-)1 



dz = -^Xi/3(r), 



(see Eq.(8.433) in 1231) 



K,{r')dr' 



cxp [— ?-cosh(z)] dz, 

cosh(z) 



(see Eq. (8.432) in [231). numerical simulations, there- 
fore, the MacDonald functions may be calculated by inte- 
grating their representations using the Simpson method, 
unless the argument, z, is very small or very large, in 
which case one can employ the series and asymptotic ex- 
pansions for cylindrical functions (see [27|). 



Appendix B. Event generator for emission 

To solve Eq. p^ numerically, one needs to pre- 
calculate the table: 

Wrnl = W {z,Xl) dz 

'x* Ix 

for discrete xi ^-i^d for discrete uniformly spaced values 
of {u}' /£)m < 1- For known xi > X* the positive value of 
Wmi most close to rnd/(5 determines the value of {ui' /£)m- 

For the ab sorp tion probability we use an analogous 
formule from l25i: 



dW. 



KKy,{r^) - J K,/,{r')dr' 

rx ^ 



£ \ d{ujt) 



UJ' 





■3 


^^3 




j z sin 


-r 
2 


V3 " 


-)1 



+00 



[ 1 • 


■3 


^^3 




/ — sm 


-r 




-)] 


J z 


2 


V3 " 





dz = --^A'2/3(r), 



dz 



2 

V3 



+ 00 



Ki,^{r')dr\ 



where 



X-it 



Here, £ is the en- 



ergy of the electron created in the 7-photon absoption to- 
gether with a positron of energy uj' — £. The x-parameter 
for a photon, x-ii is calculated as for an electron with an 
energy, w', and a dimensionless velocity, n, in terms of 
the local electromagnetic field intensities. 
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